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We generalize the level set approach to model epitaxial growth to include thermal detachment 
of atoms from island edges. This means that islands do not always grow and island dissociation 
can occur. We make no assumptions about a critical nucleus. Excellent quantitative agreement is 
obtained with kinetic Monte Carlo simulations for island densities and island size distributions in 
the submonolayer regime. 



I. INTRODUCTION 

At the present time, there is no practical approach 
to epitaxial growth modelling that bridges the gap be- 
tween microscopic and macroscopic length scales. Rate 
equations offer some hope,Ej but their reliance on uncon- 
trolled mean-field approximations remains a serious ob- 
stacle. Kinetic Monte Carlo (KMC) simulations are very 
popularp but scale-up to the micron range is very doubt- 
ful, even with future supercomputers. 

One approach to the multiple-scale problem uses an 
atomic description in the vertical (growth) direction and 
a continuum description in the lateral directions. Specif- 
ically, the random walk of individual atoms on a flat ter- 
race is replaced by the solution of a diffusion equation 
for the njipnomer density on each terrace. This is not a 
new idealj, but its recent jcebirth in the context of the 
level set (LVST) methodou is particularly promising in 
light of the relatively low computational cost needed to 
treat arbitrarily complicated surface morphologies. So 
far, good success has been achieved for irreversible epi- 
taxial growth where LVST calculations quantitatively re- 
produce the results of KMC calculations for the distri- 
bution [Of two-dimensional islands in the submonolayer 
regime. lI 

The purpose of this paper is to extend the LVST 
method to the case of reversible epitaxial growth where 
thermal detachment of atoms from island edges is al- 
lowed. This step is necessary if one hopes to produce 
a model that is relevant to growth at elevated tempera- 
tures. Moreover, a reversible LVST growth model has sig- 
nificant computational advantages over a reversible KMC 
model. This is so because KMC keeps track of every de- 
taching atom, including those that eventually return to 
the island from whence tljf y came. Such events leave the 
system unchanged overalH and slow down the simulation 
significantly. By contrast, the reversible LVST scheme we 
develop below replaces these events by their time aver- 
age and so includes only those detachments that do not 
lead to subsequent reattachments. Moreover, because 
of the mean field approach, a large number of detach- 
ment events can be treated within a single simulation 
time step. 



II. METHOD 

A. Level Sets 

The level set methodi models the time evolution of ar- 
bitrarily shaped objects in n dimensions that can undergo 
topological changes. In this paper, the relevant objects 
are two-dimensional islands and topological changes oc- 
cur due to nucleation, dissociation, and coalescence. The 
key idea is to represent a curve or interface F in 3?" by 
the level /c of a function 0(x, t) 

Tk = {x : 0(x,t) = k} x e 5R" . (1) 

Here, F^ is the set of closed curves that constitute the 
perimeters of the islands with height {k + l)a {a is the 
lattice parameter). 

Fig. |l| illustrates the level set description of a typical 
epitaxial growth scenario. The left panel is a side view 
of two islands on a terrace (a) that grow to a precoales- 
ence state (b) and subsequently merge (c). Later, a new 
island nucleates on top (d). The right panel shows the 
corresponding level set functions (p. Note that it is not 
(j) that represents the surface morphology, but only the 
level sets (cj) = and 0=1). 

The motion of F is partly deterministic and partly 
stochastic. There is a deterministic part because a mean 
field theory is sufficient to model the time-average of 
many of the physical processes that contribute to growth. 
These effects are built into a velocity function w„(x, i) 
that evolves the function (/)(x, t) in time according to the 
partial differential equation 

(/.(x,i) + i;„(x,<)|V0(x,t)| =0. (2) 

As the notation suggests, Vn is the component of the 
growth velocity in the direction of the local surface nor- 
mal n = V0/|V0|. The stochastic motion of F is associ- 
ated with nucleation events and small-island dynamics. 
There is no unique algorithm to incorporate these effects 
into (j){x,t). The particular choice we make is explained 
in detail below. 
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FIG. 1. Schematic illustration on mapping island config- 
urations during growth (left panel) onto a LS function (p. 

B. Deterministic Evolution 

It is convenient to write the velocity u„ in Eq. (|^) in 
the form 



,,dct 



(3) 



where v^^^ accounts for attachment processes that grow 
islands and w'^'^' accounts for detachment processes that 
shrink islands. The first of these is proportional to the 
diffusive flux of atoms that approach an island edge from 
its bounding terraces. Therefore, if D is the surface dif- 
fusion constant and p(x, t) is the adatom density, mass 
conservation gives 



^att ^ „2 ^ (|P 

on 



terrace 



I top of island 



(4) 



We compute the required density from the mean-field, 
driven, diffusion equation 



p(x,i) = DV2p(x, t) -f F - 2 



dt 



The loss term in Eq. 
dt 



Dai / p{y^,tfd^ 



(5) 



(6) 



accounts for the dimers that nucleate as a result of bi- 
nary collision between monomers. The total simulation 
area is aud a\ is the so-called "capture number" for 
an adatom.t2l We solve Eq. (|^) subject to the boundary 
condition that p = at every point on F. This differs 
from the boundary condition usually used for reversible 
aggregationtll because we have elected to incorporate all 



detachment effects into the velocity u'^*'* (see Appendix 
A). 

To find an explicit expression for w'^''' , we note first that 
most particles that detach from an isJand are driven back 
to that island by the diffusion field.O Our interest here is 
those particles that detach without subsequent reattach- 
ment, i.e., those that escape from the "capture zone" of 
the island. This is so because, by definition, the adatoms 
in the capture zone of a given island are guaranteed to 
attach to that island eventually. A relevant quantity is 
thus Pose J the probability that a detached particle reaches 
the border of the capture zone. The positions of these 
borders are easily calculated (the locus of points where 
Vp = 0) and so is Pcsc- We find (Appendix B) that 



log((i?i3 -f a)/flis) 

l0g(i?cz/^is) 



(7) 



where i?is and i?cz are the radii of the circularly averaged 
island and capture zone. 

We now define an effective escape rate per unit length 
of island perimeter as 



-Rdot = -Ddct Pcsc A 



(8) 



where Ddct is an effective detachment rate and A is the 
linear density of detaching particles (singly coordinated 
edge atornsl. We use A = 2 for "small" islands and the 
expressionllj 
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for "large" islands. The distinction between "large" and 
"small" islands wiUJpe made clear below. Dedg is the edge 
diffusion constant .113 From i?dct we get the desired expres- 
sion for the detachment velocity that enters Eq. (^): 



,dot 



a i?dct — a Ddct Pcsc A 



(10) 



Significantly (see below), there is negligible extra com- 
putational overhead needed to incorporate detachment 
in this way. 

It remains only to find a home for the atoms that es- 
cape from the islands. In the spirit of the mean-field 
approximation, we return them to the adatom pool by 
simply augmenting the external flux. That is, the vari- 
able F in Eq. (|) is 



F = Fo + F,, 



where _Fo is the deposition flux and 



(11) 



(12) 



is the escape rate of atoms from all island edges. In this 
integral, F runs over all level sets of (/)(x, t) (see Eq. (0)). 
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C. stochastic Evolution 

A nucleation event occurs when the variable N^uc hi 
Eq. (H) becomes larger than the next integer. This im- 
plies that (/)(x, increases by a discrete amount at a 
discrete point. In the interest of numerical stabilty, we 
smooth out this increase over several points on the nu- 
merical gid used to solve Eq. (^). The exact position 
where the dimer nucleates is chosen randomly with the 
integrand of Eq. (||) as a weight factor E 

Randomness is also important for detachment from 
"small" islands that consist of only a few atoms. Our 
approach is to choose an island area ^cut and treat all 
islands smaller than this size statistically. Thus, in a 
given time interval, we use Eq. ^ to calculate the total 
number of adatoms that detach from all islands smaller 
than Acut- This corresponds to a total area-loss Ajoss- 
If ^loss > Ai, the area occupied by one atom, we de- 
tach an atom from one of the islands smaller than ^cut- 
The specific island that loses an atom is chosen randomly 
with Pose as a weight factor. We then decrement ^loss by 
Ai and repeat the process until Ajoss < Ai. This value 
is stored and added to the loss that occurs in the next 
time interval. If a detachment process leads to an is- 
land smaller than a dimer, we dissociate the dimer and 
decrement A\oss accordingly. 

We emphasize that A^ut is not the area of a critical 
nucleus, i.e., an island that is absolutely stable against 
breakup. Instead, Acut is merely the size of the smallest 
"large" island which we use as a parameter to switch 
between statistical and continuous detachment. In the 
context of our approach, the critical nucleus is defined 



by the condition 



,,dGt 



(see Eq. (i). 



III. RESULTS 
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FIG. 2. Influence of Acut on the adatom density p (a) 
and the island density A'' (b) as a function of coverage for 
Ddet/D = 0.001. Data have been averaged over 10 runs. 



Indeed, as islands grow bigger, fewer particles escape 
from the island boundaries since the number of escaping 
particles relative to the perimeter of the growing islands 
becomes smaller and smaller. 



A. Parameter Choices & Systematics 



All the LVST simulations reported here use the value 
D/Fq — 10^. Dcdg in Eq. (^) was chosen as 10^ and 
calculations where carried out on a 200 x 200 lattice rep- 
resented on a numerical grid of 568x568 points. To de- 
termine Acut, we compared runs with different choices for 
this parameter and looked for stabilization of the physical 
results. Thus, Fig | compares island and adatom densi- 
ties as a function of coverage 6 obtained for Acut = 4,6, 
and 8 at a detachment rate of D^ct/D = 0.001. Based 
on this data and related statistical tests, we find that our 
results are independent of Acut if Acut > 6. Therefore we 
use Acut = 6 in all subsequent simulations. 

Fig. H illustrates the role of detachment during growth. 
It shows Fiov from Eq. ( pT| ) as a function of Fot for differ- 
ent values of Ddet/D. The fact that i^rev has its highest 
values at the earliest stages of growth (for the higher 
values of D^ct/D) shows that most of the effective de- 
tachment takes place when the islands are quite small. 



B. Comparison to KMC 

As a critical test, we compajed our LVST results di- 
rectly with KMC simulationsE'O In our KMC simula- 
tions adatoms are allowd to hop to a nearest-neighbor site 
at a rate r„ = D exp{—nEN / k^T) , where n is the num- 
ber of nearest-neighbours, /cb the Boltzmann-constant, 
and T the temperature. Adatoms are deposited at a rate 
_Fb. In all simulations presented in this paper the ratio 
D/Fq is set to 10^. The energy barrier i?jv is chosen 
such that Ddot — D exp{—Ei^ /kBT). In addition, singly 
coordinated edge-atoms are allowed to diffuse along the 
step edge at a rate Dcdg- We chose Dcdg=Ddct, but the 
results were not sensitive to this parameter as long as the 
islands were compact (as our LVST model assumes). 

Fig. ^, Fig. H, and Fig. ^ respectively show the adatom 
density, island density, and island size distributon ob- 
tained by both methods as a function of detachment 
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FIG. 3. Escape rate of atoms from all island edges F^cv 
(see Eq. (|l^) as a function of coverage O for different values of 
Ddct/D- Simulation parameters as in Fig. ^, except Acut = 6. 



rate. These-carsies agree with previous reversible KMC 
simulationalJllj'Ea that discuss, e.g., the physical origin 
of the observed saturation of the island densities and the 
sharpening of the island size distributions. Evidently, 
there is semi-quantitative agreement between LVST and 
KMC. The size distributions results are particularly no- 
table because they reflect information about spatial cor- 
relations that are averaged over to get p and N. 

The only disagreement we find between LVST and 
KMC is for the saturation value of the island density 
for the smallest value of the detachment rate. We under- 
stand this based on recent researchtD with irreversible 
growth where a corresponding disagreement vanishes 
when the p = boundary condition for Eq. (|^) is applied 
not at the true perimeter of each island (as we do) but 
instead at a closed boundary that exceeds the perimeter 
everywhere by one lattice constant. The disagreement 
disappears altogether when the island density is small. 
This is consistent with our observations because we get 
good agreement with KMC when the detachment rate is 
large (and the island density is small). 

Fig. 1^ shows some systematic features of reversible 
growth as a function of the effective detachment rate. 
The LVST data (collected at 0.25 ML coverage) show 
that the average adatom density increases linearly with 
_DdGt while the average island density decreases exponen- 
tially with Ddct- Unfortunately, we have been unable to 
derive theses interesting results analytically using rate 
equations. 

The data in Fig. ^ was easy to obtain because, be- 
ginning with an irreversible growth simulation, the ex- 
tra computational cost to include detachment is very 
small for LVST compared to KMC. This is so because 
LVST precisely suppresses the time-consuming detach- 
ment/attachment fluctuations that occupy a KMC sim- 
ulation. Moreover, the LVST has essentially no restric- 
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FIG. 4. Adatom densities as obtained by the level set 
method (a) and KMC (b). In the LVST calculation all pa- 
rameters are as in Fig. H, except Acut ~ 6. 



tion on the number of detachment events simulated dur- 
ing a speciflc time step. Quantitatively, Fig. || shows 
that the LVST method requires only negligibly more run 
time to include detachment whereas the KMC simulation 
cost increases sharply as the detachment rate increases. 
The LVST results depend very weakly on the rate of de- 
tachment because the increased cost is associated wholly 
with reductions in the step advance rate. Speciflcally, the 
adatom density rises with increasing detachment rate so 
the gradients evaluated in Eq. (|^) become larger. But 
the step advance rate is limited by the condition that the 
boundary of the level set function can advance only one 
grid point in each simulation step. This implies that the 
scaling should be even better for simulations of, say, an- 
nealing processes, where there is no deposition flux and 
the adatom density is very low. 
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FIG. 5. Island densities N as obtained by the level set 
method (a) and KMC (b). Parameters as in Fig. ^ 



IV. CONCLUSION 

In summary, we have developed a method to model 
epitaxial growth including atomic detachment from is- 
land edges within the context of the level set method. 
By all reasonable measures, the results are in excellent 
agreement with KMC simulations. Moreover, the LVST 
simulations scale significantly in CPU-time demand than 
KMC simulations when the effective detachment rate is 
large. This is so because our mean field method elim- 
inates the many atomic detachment events (each pro- 
cessed separately in KMC) that do not lead to successful 
escape from an island. 
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FIG. 6. Island size distributions where Us is the density 
of islands of size s, (s) is the average island size, and O is 
the coverage. Closed circles: Level set result, open squares: 
KMC. Detachment rates are in Fig. (a): Ddct/i3=0.0001, in 
Fig. (b): L»det/£'=0.0005, and in Fig. (c): L>det/O=0.001. 
Data have been sampled at O = 0.25. Other parameters as 
in Fig. |. 



APPENDIX A: BOUNDARY CONDITIONS 

The usual boundary condition for reversible growthS 
sets p = pcq at every island perimeter. In this Appendix, 
we relate poq to the detachment velocity v'^°*^ used in this 
paper. The idea is to consider two adatom densities: p 
and p as schematically shown in Fig. ||. p (p) is the sys- 
tem where the adatom density drops to pcq (p=0) at the 
island boundary. We then write down diffusion for both 
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FIG. 7. Equilibrium adatom density p (a) and equilibrium 
density of islands A*' (b). Data have been sampled at 0.25 
coverage. Other parameters as in Fig. ^ 

adatom densities and derive from those the respective 
reversible growth velocities. 

Let p (the solid line in Fig. W) be the exact solution of 



= L'V^pCx) + Fa X G 
p(x) = pcq X e dB{t) 



(Al) 
(A2) 



where D, is the domain and dBit) the island boundaries. 
Pcq is the (finite) equilibrium value of /o(x) at the island 
boundaries due to detachment of atoms. 

If C{t) is the capture zone and dC{t) its boundary 
(dashed line in Fig. |^), then with the boundary condition 



dp 
dn 



= 



x e dC{t) 



(A3) 



we can uniquely solve the diffusion equation (Al). 

In this case the reversible growth velocity is given by 
Eq. (|) (labeled Vi in Fig. |). 

Now we want to replace p by a adatom density p 
(dashed line in Fig. ^ such that 

= DV^p{x) + Fo + Frev X e O (A4) 
p(x) = X G dB{t) (A5) 
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FIG. 8. Scaling of the level set and KMC method as 
a function of detachment rate. Runtimes t are normalized 
to the runtime of the irreversible case to. Parameters as in 
Fig. |. 



were F^ev is as explained in Eq. (|11| ) the difi^usive flux 
away from the isla nd bound ary due to detachment. 
It follows from (^) and (M) that 



V^p(x 
and therefore 



- ^° + ^'-^^V^p(x) 



Fn 



P(x) = ^^^^-^r^(p(x) - Pcq) 



(A6) 



(A7) 



If we want the two systems p and p to be equivalent, then 
we must require 



p(x) (Fx — / p{x) dPx 

C(t) Jc(t) 



(A8) 



Combining (A?) and (A8) we obtain 
Fo + F,_ 

Jc{t) 



I [p{x)-p,^]d'x^ f -p[x)Sx (A9) 
^0 Jc{t) Jc{t) 



which is equivalent to 



Frev / -p(x)Sx = {Fo + F,cv) / Pcqrf'x . (AlO) 

Jcit) Jc(t) 

If A{C) is the area of the capture zone, the relation be- 
tween Frev and peq is 



1 f Pi^M _ 1 

A(C) JC(t) Pcq " ^ 



(All) 



Assuming that -Rdct is independent of dV we can rewrite 
Eq. (|l|) as 
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^rcv — -^RdctL 

Then, Eq. @ reduces to 



(A12) 



,,dct 



2* 



L 



a 



Fo 



A(C) Jc(t) p„, " 



(A13) 



APPENDIX B: DERIVATION OF Peso 

Any atom that detaches from an island but does not 
reach the capture zone boundary is driven by the diffu- 
sion field back to the island. In such a case, pcsc is zero. 
But if the detached adatom does reach the border of the 
capture zone, Pf.sc is one. Therefore, to determine the 
probability distribution P(x) for all possible paths of an 
detaching adatom, we solve the diffusion equation 



V2p(x) = 



with the boundary condition 



P(x) = 



X G 5c. 



(Bl) 



(B2) 



where Sis and S'cz are the border of the island and of the 
capture zone respectively. These boundary conditions 
represent the either successfully escape (probability 1 to 
reach S'cz) or the reattachment to the island (0 at Sis). 



Eq. (Bl) and (B2) can in principle be solved for any is- 
land and capture zone geometry. But in our calculations 
we assume for simplicity a spherical average for both, the 
islands and the capture zones. For this case the general 
solution to Eq. (il|,il) is 



F(x) = 



log(|x|/Pis) 

log(i?cz/i?i.) 



(B3) 



where R{s and i?cz are the radii of the island and of the 
capture zone. From this we obtain the escape probability 
of a detached atom by taking |x| = R{s -I- a (a is the lat- 
tice parameter) because when an adatom detaches, it is 
roughly at distance a from the island boundary. There- 
fore 



log((Pis + a)/Pis) 

l0g(i?cz/i?is) 



(B4) 
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FIG. 9. Schematic illustration to the derivation of the 
equivalence of w'^^' and peq- 
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